4. Inversione delle armoniche pure

filtri
FFT
Autore/Autrice

Mariolino De Cecco e Paolo Bosetti

Data di Pubblicazione

15 maggio 2025

Sommario

Questo esercizio genera delle armoniche pure in ingresso ad un sistema di isolamento delle vibrazioni, simulazione dell’uscita e—mediante inversione della dinamica—stimare quella originaria in ingresso.

1 Parte fornita per l’esercitazione

1.1 Funzioni

# FUNZIONE: Bodeplot
ggbodeplot <- function(tf, fmin=1, fmax=1e4, df=0.01) {
  # vector of points for each order of magnitude (OOM):
  pts <- 10^seq(0, 1, df) %>% tail(-1)
  # vector of OOMs:
  ooms <- 10^(floor(log10(fmin)):ceiling(log10(fmax)-1))
  # combine pts and ooms:
  freqs <- as.vector(pts %o% ooms)
  # warning: bode wants pulsation!
  bode(tf, freqs*2*pi) %>% {
    tibble(f=.$w/(2*pi), `magnitude (dB)`=.$mag, `phase (deg)`=.$phase)} %>%
    pivot_longer(-f) %>% 
    ggplot(aes(x=f, y=value)) +
    geom_line() +
    scale_x_log10(
      minor_breaks=scales::minor_breaks_n(10), 
      labels= ~ latex2exp::TeX(paste0("$10^{", log10(.), "}$"))) +
    facet_wrap(~name, nrow=2, scales="free") +
    labs(x="frequency (Hz)")
}

# FUNZIONE: Generazione segnale sinusoidale con diverse armoniche
signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*cos(t*f[i] + phi[i] ))))
  })
}

1.2 Segnale in ingresso

Definiamo un segnale d’ingresso (nel caso di sistemi di misura si tratta di un misurando) costituito da segnali armonici, più un eventuale disturbo normale:

f0 <- 10

# Provare con un'armonica pura
# pars <- tibble(
#   w = c(1),
#   f = c(f0),  # 1 funziona
#   phi = c(0)
# )

# Provare con tante armoniche
pars <- tibble(
  w = c(1, 0.1, 0.3),
  f = c(15, 20, 45),
  phi = c(0, 0, 0)
)

Nt <- 1000 # numero totale dei campioni
s <- tibble(
  t = 0:Nt * 0.001, # 1 kHz di frequenza di campionamento
  u = signal(t, pars), 
  un = u + rnorm(length(t), 0, pars$w[1] / 10)
)

La trasformata di Fourier mostra i picchi attesi:

s %>% 
  mutate(
    f = 0:(n()-1)/max(t),
    fft_u = fft(un),
    intensity_u = Mod(fft_u) / n()*2,
    phase_u = Arg(fft_u)/pi*180
  ) %>% 
  slice_head(n=as.integer(nrow(.)/2)) %>%
  ggplot(aes(x=f, y=intensity_u)) +
  geom_spoke(aes(y=0, radius=intensity_u, angle=pi/2)) +
  geom_point()

FFT del segnale armonico di riferimento

1.3 Impieghiamo un sistema per isolamento da vibrazioni per simulare l’uscita

Con il metodo delle impedenze generalizzate si ottiene:

\[ H(i\omega)=\frac{V_\mathrm{out}(i\omega)}{V_\mathrm{in}(i\omega)}=\frac{C i\omega + K}{M (i\omega)^2 + C i\omega + K} \tag{1}\]

La frequenza naturale del sistema è \(f_0=\frac{1}{2\pi}\sqrt{\frac{K}{M}}\), e l’attenuazione comincia a \(\sqrt{2}f_0\).

Possiamo definire la funzione di trasferimento in Equazione 1 con la funzione control::tf(), che prende come due argomenti due vettori con i coefficienti della Equazione 1, in ordine decrescente di grado della variabile \(i\omega\):

M <- 10
K <- 1000
C <- 50

# Frequenza naturale:
fn <- 1/(2*pi) * sqrt(K/M)

num <- c(C, K)
den <- c(M, C, K)

H <- tf(num, den)

ggbodeplot(H, fmin=0.1, fmax=100) +
  geom_vline(xintercept=c(1, sqrt(2)) * fn, color="red", linetype=2) +
  labs(title=paste(
    "Natural frequency:", round(fn, 2), "Hz",
    " - Isolation: >", round(sqrt(2)*fn, 2), "Hz"))
Registered S3 methods overwritten by 'signal':
  method         from   
  print.freqs    gsignal
  print.freqz    gsignal
  print.grpdelay gsignal
  plot.grpdelay  gsignal
  print.impz     gsignal
  print.specgram gsignal
  plot.specgram  gsignal

Bode plot per il sistema di isolamento da vibrazioni
cat("Per la frequenza", f0, "Hz abbiamo: \nModulo:", Mod(freqresp(H, 2*pi*f0)), "\nFase:", Arg(freqresp(H, 2*pi*f0))*180/pi, "°\n")
Per la frequenza 10 Hz abbiamo: 
Modulo: 0.08539786 
Fase: -102.9892 °
Versione alternativa

O anche, usando glue, è più semplice comporre delle stringe inserendo tra graffe le espressioni da valutare:

glue("Per la frequenza {f0} Hz abbiamo: \nModulo: {Mod(freqresp(H, 2*pi*f0)) %>% round(3)}\nFase: {(Arg(freqresp(H, 2*pi*f0))*180/pi) %>% round(3)}°")
Per la frequenza 10 Hz abbiamo: 
Modulo: 0.085
Fase: -102.989°

Dunque, se inseriamo un’armonica a fase nulla come un coseno \(\cos(\omega_{0} t)\) di frequenza pari a 10 Hz avremo in uscita la stessa armonica attenuata di 0.085 e sfasata di 103°. Dovrebbe comportarsi effettivamente come un sistema di isolamento delle vibrazioni.

1.4 Simulazione uscita

Andiamo a verificarlo. Simuliamo cosa succede in uscita se imponiamo il segnale al sistema d’isolamento impiegando il simulatore integrato in R lsim().

output <- lsim(H, s$un, s$t)
s <- s %>% 
  mutate(
    # Simulazione dell'uscita
    y = output$y[1,]
  )

# Grafici
pp <- plot_ly() %>%
  add_lines(s$t, s$y, name = "output", line= list(color= "red")) %>%
  add_lines(s$t, s$un, name = "input", line= list(color= "blue")) %>%
  layout(title = "Input & Output")
pp
Versione alternativa

Attenzione: i grafici che usano plotly sono oggetti Javascript e in quanto tali non sono compatibili con LaTeX, quindi non è possibile generare anche la versione pdf se il documento contiene grafici plotly.

(s %>% 
  select(t, input=y, output=un) %>% 
  pivot_longer(-t, names_to = "segnale", values_to = "valore") %>% 
  ggplot(aes(x=t, y=valore, color=segnale)) +
  geom_line()) %>% 
  ggplotly()

Se avete usato come parametro frequenza della funzione creata 10 Hz potete verificare come il segnale in uscita si sia effettivamente attenuato di parecchio!

Esercitazione

Riconoscimento punti esame

Il completamento di questo esercizio comporta il riconoscmiento di ±1 punto.

Da fare

Stimare il segnale in ingresso originale a partire dal segnale in uscita e dal modulo e la fase della funzione di trasferimento

NOTA: impiegare funzioni che processano automaticamente tutte le armoniche che sono presenti nel segnale. In altre parole impiegare nelle funzioni che costruirete la struttura della funzione signal()

Suggerimenti

  • calcolare FFT
  • trovare i picchi FFT
  • applicare inversa di modulo e fase della funz trasferimento alle armoniche corrispondenti ai picchi del modulo della FFT
  • scrivere una funzione che stima l’ingresso da generiche componenti armoniche in uscita che avrete salvato in un contenitore dati ‘picchi’ ed H(w): signal_input <- function(t, picchi, H) {}
  • quando il vostro codice funziona provate ad aggiungere altre armoniche
  • quando il vostro codice funziona provate ad aggiungere rumore

2 OMETTERE E FAR FARE COME ESERCIZIO

2.1 Funzioni da costruire

# FUNZIONE: calcola il modulo di una funzione di trasferimento
# w è la pulsazione in radianti
# num e den numeratore e denominatore in ordine di esponente decrescente
valore_H <- function(num, den, w) { 
  sum(map_vec(seq_along(num) , \(n) rev(num)[n]*(1i*w)^(n-1) )) / sum(map_vec(seq_along(den) , \(n) rev(den)[n]*(1i*w)^(n-1) ))  
}


# FUNZIONE: Stima ingresso da componenti armoniche in uscita 'picchi' ed H(w)
signal_input <- function(t, picchi, H) { 
    map_dbl(t, \(t) sum( map_vec(seq_len(nrow(picchi)) , \(i) picchi[i,1] / Mod(valore_H(num, den, picchi[i,2])) * cos( picchi[i,2]*t + picchi[i, 3] - Arg(valore_H(num, den, picchi[i,2])) ))))
}

2.2 Calcoliamo FFT ed i suoi picchi

# Aggiungiamo l'FFT
s <- s %>% 
  mutate(
    f = 0:(length(t)-1)/max(t),
    fft_y = fft(y - mean(y)),
    intensity_y = Mod(fft_y) / n()*2,
    phase_y = Arg(fft_y)/pi*180
  )

plot_ly() %>%
  add_lines(x = s$f, y = s$intensity_y, type = "bar", name = "FFT uscita", marker = list(color= "red"))  %>%
  layout(title = "Modulo FFT Output")
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
plot_ly() %>%
  add_lines(x = s$f, y = s$phase_y, type = "bar", name = "FFT uscita", marker = list(color= "red")) %>%
  layout(title = "Fase FFT Output [°]")
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
# Troviamo i primi N picchi del modulo della FFT
N <- length(pars$f)
picchi <- matrix(0, nrow = N, ncol = 3)  # Inizializza la matrice dei picchi
peaks <- findpeaks(s$intensity_y[1:(length(s$intensity_y) %/% 2)], nups = 1, ndowns = 1, sortstr = TRUE)

# Controlla che ci siano picchi previsti
if (!is.null(peaks) && nrow(peaks) >= N) {
  picchi_values <- peaks[1:N, 1]  # Valori dei picchi
  picchi_indices <- peaks[1:N, 2] # Indici dei picchi
  picchi_freq <- s$f[picchi_indices]  # Frequenze corrispondenti ai picchi
  picchi_fasi <- s$phase_y[picchi_indices] # Fasi corrispondenti ai picchi

  # Assegna i valori alla matrice picchi
  picchi[, 1] <- picchi_values
  picchi[, 2] <- 2 * pi * picchi_freq
  picchi[, 3] <- pi * picchi_fasi / 180
}
print(picchi)
            [,1]       [,2]       [,3]
[1,] 0.054318212  94.247780 -1.7265506
[2,] 0.006344690   6.283185  0.1546009
[3,] 0.005216083 282.743339 -1.6356877

2.3 Applichiamo inversa della funz trasferimento

# Applichiamo inversa della funz trasferimento ai picchi della FFT in modo da stimare l'ingresso a partire dall'uscita
s <- s %>% 
  mutate(
    u_est = signal_input(t, picchi, H)
  )

pp %>%
  add_lines(s$t, s$u_est, name = "input stimato", line= list(color= "green")) %>%
  layout(title = "Input originario e stimato")